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0> ■ Abstract 

We study surface effects in amorplious polymer systems by means of com- 
■ puter simulation. In the framework of molecular dynamics, we present two 

different methods to prepare such surfaces. Free surfaces are stabilized solely 
by van-der-Waals interactions whereas confined surfaces emerge in the pres- 
' ence of repelling plates. The two models are compared in various computer 

simulations. For free surfaces, we analyze the migration of end-monomers to 
the surface. The buildup of density and pressure profiles from zero to their 
bulk values depends on the surface preparation method. In the case of con- 
\^ I fined surfaces, we find density and pressure oszillations next to the repelling 

' plates. We investigate the influence of surfaces on the coordination num- 

"c^ ' ber, on the orientation of single bonds, and on polymer end-to-end vectors. 

^ . Furthermore, different statistical methods to determine location and width 

, of the surface region for systems of various chain lengths are discussed and 

', applied. We introduce a "height function" and show that this method allows 



o 



o 



to determine average surface profiles only by scanning the outermost layer of 
monomers. 
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INTRODUCTION 



The technological applicability of polymeric materials depends in a number of important 
cases on the polymer's surface properties: the surface roughness determines the friction co- 
efficient, cracks emerging from the surface influence the material's durability, and surface 
defects reduce the transparency of polymeric panes. At the same time, the experimental 
and manufacturing techniques advance to smaller and smaller scales. Nowadays one can 
pursue mechanical effects on the nano-scale. This opens new, challenging, and exciting 
roads for academic research and promising prospects for industrial applications. Nano- 
mechanics influences the design of computer chips and is fundamental for new developments 
like miniaturized optical devices and nano-engines. A related topic is the mechanics of very 
thin films with thicknesses in the nano-range. But again, it is not just the bulk behaviour 
that decides about the nano-technological applicability of a certain material but its surface 
properties. The nano-scale processes that lead to macroscopically relevant surface-related 
phenomena like friction, wear, hardness, and certain types of material failure like surface 
degradation are scarcely understood from the theoretical point of view. In some cases, 
phenomenological theories exist which attempt to systematize the various experimental 
observations. A fundamental understanding, however, tentatively based on fundamental 
nano-scale processes, is lacking. The reason is that on this scale, the number of mutually 
interacting particles is very high, but the system is still too coarse to apply continuum 
mechanics in a straightforward way. Some simplifications can be made in the case of crys- 
talline solids but amorphous systems like polymer glasses still represent a major challenge. 
Techniques from computational physics can help to deal with the many degrees of freedom 
involved 0. 

In ongoing work, we intend to address several basic questions that pertain to the un- 
derstanding of nanoscopic amorphous systems in conjunction with tribological and surface 
effects like friction, wear, and hardness. Our investigation is mainly restricted to polymeric 
materials. Crystalline solids have been studied in a beautiful series of pioneering papers 
by Landman et al. [Q. In the present contribution, we concentrate on the preparation of 
polymer surfaces in the framework of computer simulation. We are furthermore concerned 
with the detailed characterization of important surface properties. This is a necessary foun- 
dation to correctly interpret our findings from a large number of dedicated polymer surface 
indentation simulations 0. These simulations provide a prototype nano-scale experiment 
to probe surface hardness, friction, and damage mechanisms. Related to these activities are 
computer studies on the deformability of nano-scale polymeric films . 

POLYMER AND SURFACE MODELLING 

For the polymer chains, we use a united atom model in conjunction with Newtonian 
dynamics. In addition to harmonic chain forces which keep the bond lengths next to the 
equilibrium value, we model the fiuctuation of bond angles, again by a quadratic potential. 
Between monomers which do not participate in mutual bond length or bond angle interac- 
tions, Lennard- Jones forces are acting, both to model an excluded volume effect and to hold 
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TABLE I. Parameters of the polyethylene model. 


Lennard-Jones energy, e 


8.3027 • 10^22 J 


5.18 meV 


Lennard-Jones length, a 


380 pm 


3.8 A 


monomer mass (CH2 group), m 


2.3248 • 10 ^ kg 


14 atomic units 


unit of temperature, e/Zcg 


OU.ioO/ IV 




LllilL Ui liidoo U.cllDlL_y, IlilU 


^ZiO.UUOl iS.^/ Ill 


U.^^O/ ^/ Clll 


unit 01 time, [ma /e) 


2.0108 ps 




., f» -| ., / / \]_/2 

unit of velocity, [e/m) ' 


188.9822 m/s 




unit 01 rorce, e/a 


2.1849 pN 


1.36 meV/A 
0.36 meV/A 


unit 01 spring constant, e/a 


5.7498 • 10 J\ /m 


unit of pressure, e/cr^ 


151.3103 bar 


0.09 meV/A 


f prjTnpT^ifnrp 'T"' 


301 K and 361 K 




bond length, Iq 


152 pm 


1.52 A 


bond angle, 


109.47° 


cos 00 = —1/3 


spring constant (bond length), kh 


5.7498 • 10^ N/m 


3.59 eV/A^ 


bending constant (bond angle), kg 


8.3027 • 10-19 J 


5.18 eV 


simulation time step, 


2.0108 fs 





the polymer system together. Note that we neglect any torsional potential in the present 
study. To be explicit, the Hamiltonian of the model is of the general form 



H = H, 



bondlength 



bondangle 



LJ 



-^bondlength 
-^^bondangle 



X! y (^bond — ^0)^ , 
bonds 

y (C0s6'angle " COS 6^0)^ , 

angles 
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The Lennard- Jones interaction is implemented with a cutoff of 2.5 a and appropriate poten- 
tial and force shifts are used to retain continuity. We use molecular dynamics methods to 
compute the motion of the monomers. Models of this kind are described at various places in 
the literature 0]. We intend to capture some essential features of polyethylene chains and 
appropriate model parameters are compiled in Tab. | (see also Ref. [§). 

We employ two alternative methods to prepare amorphous polymer surfaces. The con- 
fined polymer surface, Fig. |I], is realized by plates that repell the monomeric units and 
thereby restrict the polymer to a certain region in space. This gives rise to a surface model 
which is mainly characterized by its smoothness. In the simulations, we set up the united 
atom model with periodic boundary conditions in the x- and y-direction. In the z-direction, 
we put two repelling plates at a distance chosen such that the resulting bulk density of the 
polymer system settles to the desired value. The plates interact with the polymers through 
the repelling (r-^^) part of the Lennard- Jones potential. The other alternative, Fig. ^ is the 
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This is a dummy figure ... 



FIG. 1. The confined polymer surface (at the top of the bulk). The horizontal line of the frame 
represents the repelling plate. The monomers are drawn with their Lennard-Jones radius a. In 
the lateral directions, periodic boundary conditions apply. 



This is a dummy figure ... 



FIG. 2. The free polymer surface (at the top of the bulk). The horizontal line of the frame is 
out of the interaction range. The monomers are drawn with their Lennard-Jones radius a. In the 
lateral directions, periodic boundary conditions apply. 
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FIG. 3. Various measures of the mean square displacement (MSD). In particular, gi is the MSD 
of the center monomers, §2 is the same but relative to each chain's center of mass, 53 is the MSD 
of the chains' center of mass, 54 is the MSD of the end monomers, and 55 is the same but again 
relative to each chain's center of mass. The chain length is = 40 and only every tenth sample 
point is marked by a symbol. 

free polymer surface. Here the polymer bulk is held together solely by the overall attracting 
intermolecular van-der-Waals interactions. Individual chains may evaporate or condensate. 
The emerging surfaces are rougher than in the confined case. We suggest that the confined 
model applies to the interface between the polymer and a comparatively rigid material like 
a metal. Such surfaces may also be formed on explosion pressed materials. On the other 
hand, the free surface is capable to model the interface between the polymer and a (dilute) 
medium where any mutual interactions can be neglected (in effect, it is the interface to a 
vacuum) . 

All simulations to probe surface properties are perfomed at reduced temperatures of 
T* = T/^e/k-o) = 5.0 (confined system) or 6.0 (free system). These temperatures correspond 
to 301 or 361 Kelvin which is above the melting temperature of real polyethylene. We stress 
that our simulations are not yet supposed to fully mimic real polyethylene (it has been 
demonstrated, however, that much can be achieved by carefully tuning the parameters of 
united atom models 0). The chain lengths are still too short and we treat chain-end 
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FIG. 4. Distribution of density and end-monomers for the confined surface. Shown is the num- 
ber density of monomers (soUd curve, left axis) and the relative frequency of chain-end monomers 
(bar chart, right axis) in slices parallel to the surface. The polymer system is monodisperse and 
contains 205 chains of = 40 monomers each. The horizontal dashed line indicates the mean 
value for chain-ends: 2 ends/40 monomers = 0.05. Also shown is the location of the confining 
plates at both sides. 
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monomers like mid-chain monomers. We also neglect the torsional potential (the rotation 
about C-C axes is not restricted). The high temperature helps to ensure that in feasible 
simulation times, the chains travel several radii of gyration in order to guarantee correct 
thermodynamic properties. To check this, various mean square displacement measures as 
defined in Ref. are shown in Fig. ^ Indication that the chain in its entirety has moved 
through the matrix comes from the relative displacements of the inner monomers relative to 
the center of mass motion ((72) and the outer monomers relative to the center of mass motion 
{g5). The cross-over points define charateristic time-scales of the system. The presence of 
the cross-over points in Fig. |^ indicates that the simulation time scale is sufficiently long 
compared to the short-order relaxation times. Note that the data for all the figures to 
follow represent averages over a great number of system configurations which have been 
sampled from simulations which extend over a span of time comparable to that in Fig. |^. 
For each chain length (A^ = 20, 40, 60, and 80), the simulation follows 8200 monomers over 
1.4 nanoseconds. If not otherwise indicated, the results given are for the chain length iV = 40. 



CONFINED POLYMER SURFACES 

In the case of the polymer system with confined surfaces, the plates are rapidly brought 
from infinity to a distance such that the monomer number density in the resulting box is 
Pbuik = 2.0. The shock results in a distribution of the density in the z-direction (system 
height) as shown in Fig. ^. These typical periodic changes of the density distribution close 
to the surface |TD[ are due to the external pressure produced by the plates, which leads to 



a local partial crystallization. The same behaviour has been found in simulations analyzing 
surface phenomena and studied in the context of one long chain between two parallel 



plates |]T2[. The chains experience strong packing constraints as well as a loss in entropy 
that must show up also in the surface properties that can be observed during indentation. 
Fig. ^ also probes the chain-end enrichment effect for confined surfaces. Such an enrich- 



ment has been noticed by Wang and Binder in a lattice model. Shown with the density 



profile is the relative frequency of end-monomers. This is computed by defining slices of 
height d in the 2;-direction and dividing the number of chain ends by the total number of 
monomers within each slice. The resulting relative frequency is fairly independent of the slice 
height d. We conclude from the figure that there is no appreciable chain-end enrichment at 
confined surfaces. This is in accordance with the fact that there are no energetic or entropic 
reasons for the end-monomers to gather in the vicinity of the plates (see the discussion for 
free surfaces below). The constraints imposed by the plates can favour a high positive virial 
contribution to the pressure. Remember that for confined surfaces, the volume and therefore 
the mean density of the simulation system is prescribed. The high pressure in Fig. |^ can be 
traced back to repelling intermolecular interactions caused by dense packing. 



FREE POLYMER SURFACES 



To study the effects of a rougher and presumably more realistic surface, we set up a 
free system with periodic boundary conditions in x- and y-direction and free boundary 
conditions in ^-direction. The polymer bulk is held together by the attractive part of the 
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FIG. 5. The pressure distribution vs. the height direction for the confined surface model (pres- 
sure in reduced units). The total pressure is the sum of the virial and the kinetic contribution. Also 
shown is the location of the confining plates at both sides. The dashed horizontal line indicates 
the expected value for the kinetic (ideal) part: pf(jeai ~ P*'^* ~ ^-^ ' ^-^ ~ 10.0. 
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FIG. 6. The monomer number density in the bulk (open circles) for different chain lengths 

(20, 40, 60, and 80 monomers) plotted vs. the reciprocal chain length in the free surface 
model. The dashed line is a linear fit to the data points. Its intercept with the y-axis yields 
Pbuik(^ ^ do) « 2.12 (slope « -6.87). 
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inter-molecular Lennard- Jones interaction. At the beginning, this system can be prepared 
like a confined one with two repelling plates in the 2;-direction. This is mainly done to bring 
the polymer bulk initially into a regular shape. By previous experience, the density during 
the constrained phase is chosen close to the bulk density of a free system. After a time, the 
plates are removed and the system is free to relax into the surrounding vacuum. Another 
possibility is to start with a periodic box completely filled with polymer. This system is run 
over a relaxation phase which is finished when the chains have diffused over at least two radii 
of gyration. After that the periodic boundary condition in the z-direction is removed and the 
box is rendered three times as high as before. The polymer film is placed into the middle of 
the box and the simulation continues. The density profiles quickly approach their stationary 
shape and we made sure that their bulk or plateau values are independent from the initial 
density in the periodic box. Especially for longer chains this gives a good criterion that 
the dynamics have approached the equilibrium stage. The bulk densities found for varying 
chain lengths are summerized in Fig. ^. The reason for the different bulk densities is the 
ratio of chain-end to mid-chain monomers. In our model, the effective volume of chain-end 
monomers is higher than that of mid-chain monomers, since end-monomers with only one 
bond need more room than those with two bonds (the Lennard- Jones range a is the same 
for all monomers and longer than the bond length Iq) and also since the end-monomers 
have a higher mobility (entropic effect). For different chain lengths, the bulk densities are 
interpolated by the formula 

'^chain ''^monomer ^ ^monomer 

Khain 2V^nd + (A^ — 2)Vmid Vmij 

Since Vend — Knid > 0, there is a 1/A^ chain-end effect with a negative slope, in accordance 
with Fig. ^ Moreover, note that in the period of time accessible to our simulations, no 
chains escape the polymer bulk. Therefore the effective pressure in the vacuum phase is 
zero (see also below. Fig. p^T]) . 
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Density profile, cliain— end enrichment, and pressure 

In Fig. 1^, the distribution of the monomer number density vs. the height {z) direction 
(the direction normal to the surfaces) is plotted. As opposed to confined surfaces, no partial 
crystallization can be observed at free surfaces. Also shown in this figure is the chain-end 
enrichment. Again in contrast to confined surfaces, we observe the clear tendency for the 
chain-ends to escape from the bulk. Chain-end monomers require a greater static volume 
than mid-chain monomers. To reach close packing, they stick to the surface. This is also 
favoured by reasons of entropy, since the end-monomers possess a higher mobility (they are 
only constrained by one chemical bond). These arguments, however, do not apply if the 
surface is pressed by repelling plates. In Fig. |^ one clearly sees the enrichment effect, but 
an exact quantitative evaluation is out of reach, since if we increase the number of bins next 
to the surface, the bin width becomes too small to achieve good statistics. The presence of 
end-monomers of high mobility and the local brush-like structure caused by the chain ends 
will certainly affect the mechanics of surface disturbances like in an indentation process. 
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FIG. 7. Density distribution and chain-end enrichment for the free surface. Shown is the num- 
ber density of monomers (soUd curve, left axis) and the relative frequency of chain-end monomers 
(bar chart, right axis) in slices parallel to the surface. The polymer system is monodisperse with 
a chain length of iV = 40. The horizontal dashed line indicates the mean value for chain-ends: 
2 ends/40 monomers = 0.05. 
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FIG. 8. The relative frequency of the normalized monomer index at free surfaces for different 
chain lengths. If all monomers in a chain would appear with the same probability at the surface, 
the relative frequency would be 1, independent of the normalized index (dashed horizontal line). 
However, higher normalized indices are more frequent, a signature of the chain-end enrichment at 
free surfaces. 
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A complementary way to analyze chain-end enrichment is to look at the distribution of 
the normalized monomer index T which is defined as 



N/2 



N/2 



2i 
N 



(3) 



Mid-chain monomers get the normalized index 0, for end-monomers it is 1. Let n(T) be 
the absolute frequency of the normalized index T found at the surface and ritotai = J2i''T'{T) 
the total number of all monomers at the surface. These quantities can be determined by 
dividing the surface into patches, i.e. dividing the whole film into columns. Within each 
column, we look for the extremal monomers (in a sense, these monomers then define the 
surface). The same approach is chosen below to sample the surface height function. The 
relative, chain length independent frequency of T is then computed according to 

m = ^ = ^ f4) 

^ ' ntotJ{N/2) 2r2totai ' ^ ^ 

The result for chain lengths = 20, 40, 60, and 80 is shown in Fig. |^. If each monomer 
in the chain would appear with the same probability at the surface, the relative frequency 
would be 1, independent of the normalized index (dashed line in the figure). We observe, 
however, a significant increase in frequency starting with a normalized index of around 0.8. 
Moreover, the shape of the curves is very similar for the different chain lengths. 

If the chain-ends migrate to the surface, one might suspect some correlation between 
the orientation of the surface and the chain end-to-end vectors. This is analyzed in Fig. |^. 
For a given chain, let Ri denote the monomer positions. Then R = R^ — Ri is the end-to- 
end vector. As before, the surfaces are parallel to the x?/-plane. We want to look at the 
normalized projections of the end-to-end vector normal and parallel to the surface, 

p2 p2 _|_ p2 

cl(z) = Jll- and cf,{z) = ^ ^ (5) 

|^|2 11^ |^|2 ^ ^ 

To place these projection coefficients into bins along the z-direction for statistical analysis, 
we use the mid-point of the end-to-end vector, z = (-Rnz + -Ri^)/2. Fig. |^ shows that near 
the surface, the end-to-end vectors are on average parallel aligned to it. This indicates 
that both end-monomers of a chain are at the same surface. The alignment, however, 
does not extend into the plateau region of the density profile. This is an important result, 
since for thin films, orientational correlations could prevent the formation of a bulk zone 
at all. Furthermore, the average length of the end-to-end vector in the simulations is 19 
A. Compared to the film height of about 80 A, there is no separation of scales. A related 
measure for orientation is the bond order parameter, b{z), defined as 

where b with components b^, by, and bz denotes an arbitrary bond vector connecting two 



monomers. In Fig. 10, the average is over chain-end bonds only. The mean values are 



computed from a number of statistically independent configurations in bins (slices) of finite 
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FIG. 9. The orientation of the end-to-end vector in the free model. Plotted is the mean square 
value {cj_{z)) = {Rl/\R\'^) for the normal-to-surface component after normalization. The mean 
values are computed in slices parallel to the surfaces at different height. In the bulk, we expect due 
to equipartition of all directions that (i?2/l-^P)buik = 1/3 (horizontal dashed line). The vertical 
lines mark the extension of the surfaces as given by the 10/90-rule applied to the number density 
profile in Fig. 0. The average length of the end-to-end vector is found to be 19 A. This represents 
a lower bound for the bin width in the figure and severely limits the number of independent data 
points. 
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FIG. 10. The bond order parameter for the chain-end bonds in the free model. The orientation 
of the end-bonds is plotted vs. their distance from the system's center of mass (left axis). This 

involves an average over both surfaces of the film. Good statistics for the end-bonds is difficult to 
obtain. In chains of length N = 40, only 5 percent of the monomers are chain-ends. In addition, the 
number of all bonds becomes very small in the surface region, as can be seen in the bin occupation 
curve (right axis). 
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FIG. 11. The pressure distribution vs. the height direction for the free surface model (pres- 
sure in reduced units). The total pressure is the sum of the virial and the kinetic contri- 
bution. The dashed horizontal line indicates the expected value for the kinetic (ideal) part: 
^'fdeai ~ p*T* 2.0-6.0 = 12.0. The dot-dashed line is the zero pressure line. 



width parallel to the film surfaces {xy-plane) . The bond order parameter is plotted vs. the 
distance of the bins from the center of mass. In the bulk, we expect isotropic behaviour, 
yielding 6 = 0. We want to test, whether the chain-end bonds in the vicinity of the surface 
are oriented perpendicular to it, which means bl > b^^y resulting in 6 > 0. Indeed, from 
Fig. ITUI we infer that this is the case. However, sufficient statistical data for the orientation 
of the end-bonds at the surface are difficult to obtain. The reason is that we sample in a 
region in space where the density is already very low. 

Fig. O shows the pressure distribution. In effect, a free surface is the interface to a 



vacuum. This means that the external pressure imposed on the free polymer film is zero. 
Kinetic and virial contributions should therefore compensate each other. In the figure, the 
total pressure is slightly less than zero and the system will still contract somewhat during a 
final relaxation phase that has not been reached in the simulation. 
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Surface profiles and surface width 



In this section, we are concerned with possible analytical fits for surface profiles and the 
systematic determination of surface location and width. A first idea is to start from the 
(number or mass) density profile in Fig. |^ and to look for the points zio and zgo where 10 
and 90 percent of the bulk density are reached. The surface width is the distance between 
these two points, C, = \ zio — Zqq \ and as the surface base point we choose the the midpoint, 
= {zio + zqq)/2. In what follows, this method will be called the "10/90-rule" . A second 
possibility is to fit the density profile with some suitably parameterized analytic function 
and to use these parameters to actually characterize the surface. This allows to compare 
profiles which have been fitted using the same family of functions. An appropriate family is 
based on the hyperbolic tangent, 

P(^) = ^Pbuik ^1 - tanh —J^^ ' (7) 

which is essentially the Fermi function, since (1/2)(1 — tanhx) = (1 + exp 2a;)^^. The 10/90- 
rule applied to a tanh-fit yields the relations Zto = zq and = ^/ lii(0. 9/0.1) ~ 0.46^. Such 
a fit is used in Fig. |l^. Another possibility is to use the error function (the probability 
distribution for the Gaussian normal probability density) for the fit, 

Piz) = Pbuik r /(C) dC with fiz) = exp (- ^'~'f^' ] . (8) 

The connection to the 10/90-rule is given by z^q = zq and a ~ 2.56^. Note that these 
functions also appear in the following paragraph when characterizing the surface by use of 
a "height function", but then they enter in a purely statistical context. 

Because of their steepness, surface characterization based solely on the flanks of the 
density profiles is problematic. The profiles must be determined by sampling monomers. 
If we put more bins onto the flanks, the bin width soon becomes too small for sufficient 
statistics. The following method offers a complementary approach. We look onto the upper 
surface from above and divide it - in the simplest case by use of a regular grid - into a number 
of patches. Within each patch, (which represents a column in the three-dimensional film) 
we identify the monomer with the highest z-coordinate. In this way we obtain a mapping 
from the patches to a set of height values. We call this mapping the "height function". 
For the given set of height values, we determine mean value and standard deviation and 
compute from them the location and the width of the surface, respectively. To get the 
result for the lower surface, we have to pick the lowest monomer within each column. The 
method is related to the surface contour one obtains using a scanning force microscope. A 
resolution of 1 or 2 Lennard- Jones a for the resolution (distance between grid lines) is found 
to be appropriate. When we plot the relative frequency of height values found, we obtain a 
probability density as in Fig. which can be fitted with a Gaussian normal density. Then 
the corresponding distribution function is determined which represents some average surface 
profile. 

The relation between the height function profile and the monomer number density profile 
can be derived in the special case that the monomer density assumes its bulk value instantly 
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FIG. 12. The free surface profile as obtained by the statistical analysis of the height function. 
The circles mark the probability to find the highest (or lowest, respectively) monomer of a column 
at the given height. This probability density comes from the simulation and has been integrated 
numerically to determine the associated probability distribution function which is marked by the 
squares. The probability distribution represents the "surface profile" in the height function ap- 
proach. The solid curves are analytical fits to the data. For the probability, Gaussian normal 
densities with (Tl/r = 2.56 ^l/r are used. The surface profile is a properly located and scaled error 
function. The surface widths (^l and ^r) and boundaries (dashed lines) as given by the 10/90-rule 
are also indicated. 
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below the surface. In this case, the finite slopes in the average density profiles stem only 
from the fact that the surface is rough. The surface is defined by the height function h{x, y) 
over a domain in the x?/-plane with base area A. For the height coordinate z we obtain the 
probability density 

Pm{z) = ^ JJ dxdyS{h{x,y) - z) . (9) 

The associated distribution function is (O ( z ) is the unit step function) 

Pm{z) = f_jM{OdC=^jjdxdye{h{x,y)-z). (10) 

On the other hand, the mean density at height z is computed by "counting all columns filled 
with matter", 

P{z) = ^ ffdxdy Q{h{x, y)-z). (11) 



A 

Under the condition mentioned in the beginning, one therefore finds 

p{z) = phnlkPhiiz) . (12) 

Of course, this is what one would have guessed knowing that the density reaches its constant 
bulk value instantly below a surface with zero local thickness. In Fig. |I3| the thickness 
measures obtained from monomer density profiles and height function statistics are compared 
to each other for different polymer chain lengths. The fit ^ ~ A^~^ was successful. We find 
that for longer chains the height function method, based on the statistics over outermost 
monomers, certainly overestimates the surface width. 



Coordination number 

In connection with amorphous materials, the coordination number of a given particle 
is the number of all its partners involved in van-der-Waals (Lennard- Jones) interactions. 
From the simulations, this number can simply be determined by counting all the monomers 
in the interaction range (a sphere of radius 2.5 Lennard- Jones a) of a given monomer 
(monomers in the same chain do count except they take part in bond length, bond angle, 
or torsion interactions). Again we assume homogeneity in the lateral {x and y) directions, 
associate with every coordination number the coordinates of its monomer, and set up the 
distribution of the average coordination number vs. the height (z) direction. The result is 
shown in Fig. 0. A theoretical estimate for the coordination number distribution c{z) must 



be based on the monomer number density profile p{z). Geometric relations then yield the 
convolution type integral 

ciz)= r%(C)^ i^'-ir-Cr) dC. (13) 

For the monomer density profile, p(-2), we use the ansatz 
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FIG. 13. The surface width ^ vs. the reciprocal chain length [N = 20, 40, 60, and 80) in 
the free model. Slightly different values for ^ are obtained by applying the 10/90-rule either to 
number density profiles as in Fig. |^ (circles) or to probability distribution profiles based on the 
height function as in Fig. |l2| (squares). The plot suggests a scaling ~ A^~^. The asymptotic 
values resulting from linear fits (dashed lines) are ^ [N oo) = 6.8 A for the density profile and 
^ {N — > oo) = 8.1 A for the height function method. 
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Height [A] 

FIG. 14. Comparison between density profile and coordination number for the free surface 
model. The simulated density curve is the same as in Fig. |^ For the density fit, a combination of 
two hyperbolic tangent functions is used. This fit allows a theoretical estimate for the coordination 
number which is given in the figure, too. In the worst case, the estimate must be determined by a 
numerical integration, whereas the determination of the coordination numbers from the simulation 
is more time consumptive. Over all sample configurations, a double loop must be run, since for 
each monomer we have to count all its interaction partners. Note that the measured coordination 
number curve is even broader than the theoretical one. 
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p{z) = pbulk 



Q(-z) 1 + tanh 



Z + ZtL 



+ Q(z) 1 - tanh 



Z — ^tR 



tL 



tL 



(14) 



with the parameters a = 9.8 A, pbuik = 1-947, Zti^ = 31.16 A, Zt^ = 31. 16 A, ^tL = 4.32 A, 
and ^tR = 4.32 A. This results in the computed curve for the coordination number in Fig. |T^. 
The discrepancy in the plateau value derives from the fact that in the estimate we implicitly 
also count the neighbor monomers in the chain (which should be neglected, see above). As 
expected, the flanks of the coordination number distribution are broader than the flanks of 
the density profile. This means that in the interaction picture, the surface region is further 
extended than according to the local density criterion. 



CONCLUSION AND OUTLOOK 

We have presented two methods to prepare amorphous polymer surfaces in molecular 
dynamics simulations with united atom models for the polymer chains. The ideas involved 
can easily be transfered to both more detailed models (which, as one extreme, explicitly 
treat each atom in a monomer, or, like the ellipsoidal model [|14[ for BPA-PC, extend the 
united atom model by considering the geometric shape of a monomer) and more coarse- 
grained polymer models. The use of classical molecular dynamics algorithms is not manda- 
tory. Monte-Carlo and Hybrid Monte-Carlo techniques can be employed as well and the 
inclusion of stochastic Langevin forces (Brownian motion) is straightforward (and, in fact, 
implemented in our simulation program). The present work is also a further step in extend- 
ing the versatility of The Materials Explorer, a continuously growing, integrated software 



tool for the simulation, visualization, and quantitative analysis of engineering materials [|15 
The work presented here will allow us to analyze our simulations on nano-mechanical surface 
indentation processes in which we prepare a surface using one of the methods discussed 
above. After that, the impact of the indentation tool causes local disturbations and the de- 
viation of the certain observables from their equilibrium values gives us information about 
the surface resistance and friction properties. Related to the indentation processes are simu- 
lations on thin film bending 0. A tool, which is more carefully driven as in the indentation 
experiment, causes non-local, large-scale deformations in the polymer film. Again vari- 
ous observables (density distribution, orientation of the chain molecules) deviate from their 
equilibrium values and have to be analyzed. 
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